home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / FSULTRA1.ZIP / ULTRXCOD.INC < prev   
Text File  |  1992-06-18  |  9KB  |  375 lines
  1. comment ! 
  2. FSU - ULTRA    The greatest random number generator that ever was
  3.         or ever will be.  Way beyond Super-Duper.
  4.         (Just kidding, but we think its a good one.)
  5.  
  6. Authors:    Arif Zaman (arif@stat.fsu.edu) and
  7.         George Marsaglia (geo@stat.fsu.edu).
  8.  
  9. Date:        27 May 1992
  10.  
  11. Version:    1.05
  12.  
  13. Copyright:    To obtain permission to incorporate this program into
  14.         any commercial product, please contact the authors at
  15.         the e-mail address given above or at
  16.  
  17.         Department of Statistics and
  18.         Supercomputer Computations Research Institute
  19.         Florida State University
  20.         Tallahassee, FL 32306.
  21.  
  22. See Also:    README        for a brief description
  23.         ULTRA.DOC    for a detailed description
  24.  
  25. -----------------------------------------------------------------------
  26. ! 
  27. .386
  28. ; This file is included in after a header file which defines
  29. ; language dependent elements. (see ultra.doc for details)
  30. ;
  31. rinit    proc        ; INITIALIZE SEED ARRAY ========================
  32.     Enter2arg
  33. ;
  34. ; On entry:
  35. ;    eax and ebx contain congruential and shift-register seeds
  36. ;
  37. ; On exit:
  38. ;    The swbb array is set using the McGill Super-Duper generator,
  39. ;    which is a mix of a congruential and shift-register sequence.
  40. ;    Flags and counters are reset.
  41. ;
  42. ; Registers Clobbered:  AX,BX,CX,DX,SI,DI.
  43. ;
  44. ; Algorithm:
  45. ;    Starting from lsm of x[0] to the msb of x[N-1], each bit is
  46. ;    formed using the sign bit of the xor of a congruential generator
  47. ;    with seed ConX and a shift register sequence with seed ShrX.
  48. ;
  49. ; Register usage
  50. ;   eax contains congruential seed
  51. ;   ebx contains the shift-register seed
  52. ;   cx contains count for two loops
  53. ;      ch counts the outer loop (for each byte of the array x)
  54. ;      cl counts for each bit of x[i]
  55. ;   es:di point to where x[i] is
  56. ;   edx is a temporary register
  57. ;
  58.     mov ch,N*4
  59.     mov di,offset swbseed    ; do loop for x[0], ... , x[4*n] (bytes)
  60. nextbyte:
  61.       mov cl,8
  62. nextbit:
  63.     mov  edx,dword ptr 69069
  64.     mul  edx
  65.  
  66.     mov  edx,ebx
  67.     shr  edx,15
  68.     xor  ebx,edx
  69.     mov  edx,ebx
  70.     shl  edx,17
  71.     xor  ebx,edx
  72.  
  73.     mov  edx,eax
  74.     xor  edx,ebx
  75.     rcl  edx,1
  76.  
  77.     rcr  byte ptr [di],1   ; shift it into the answer
  78.     dec  cl
  79.     jnz  nextbit
  80.     inc  di             ; Store the answer in swbb[i]
  81.     dec ch
  82.     jnz nextbyte
  83. ;
  84. ; reset all counters, flags etc. and return.
  85. ;
  86.     xor bx,bx
  87.     mov swb32.c,bx
  88.     mov swb16.c,bx
  89.     mov swb8.c,bx
  90.     mov swb1.c,bx
  91.     mov flags,bl
  92.     mov congx,eax
  93.     Exit2arg
  94. rinit    endp
  95.  
  96. SWBfill    proc near    ; SUBTRACT-WITH-BORROW ========================
  97. ;
  98. ;  On Entry:
  99. ;     DS should point to the data segment.
  100. ;     BX points to the array where the results will be stored (swb??.c)
  101. ;
  102. ;  On Exit:
  103. ;     The swbseed array contains all new values computed by the
  104. ;     subtract-with-carry generator. The CARRY byte contains the
  105. ;     state of the carry flag due to the last subtraction.
  106. ;     swb??.x contians the seed xored with a congurential.
  107. ;
  108. ;  Registers Clobbered: AX,BX,CX,DX,SI,DI
  109. ;
  110. ;  Algorithm:
  111. ;     The following subtractions are performed from right to left.
  112. ;     The carry propagates as though these were two long numbers,
  113. ;     with each x[i] being a `digit' in base 2^32.
  114. ;
  115. ;        x[12] ...  x[ 0]      x[36] ...  x[13]
  116. ;       -x[36] ... -x[24]     -x[23] ... -x[ 0]
  117. ;       ---------------------------------------
  118. ;        x[36] ...  x[24]      x[23] ...  x[ 0]
  119. ;
  120. ;
  121. ;    for a bigger table, it could be done in three steps by:
  122. ;
  123. ;     x[ 12] ...  x[   0]  :  x[N-1 ]  x[N-24]  :   x[N-25] ...  x[13]
  124. ;    -x[N-1] ... -x[N-13]  :  x[N-14] -x[N-37]  :  -x[N-38] ... -x[ 0]
  125. ;    ----------------------:--------------------:---------------------
  126. ;     x[ 36] ...  x[  24]  :  x[13  ]  x[0   ]  :   x[N-1 ] ...  x[37]
  127. ;
  128. ;     The x's also could be considered as pairs of base 16 digits,
  129. ;     so that x[i] is the pair y[2i+1]y[2i]. This allows us to use
  130. ;     only 16 bit subtractions with carry, perfectly suited for all
  131. ;     80x86 processors. The same idea could be extended for machines
  132. ;     with only eight bit, or even only 1 bit arithmetic.
  133. ;
  134.     EnterFill
  135.     mov ah,byte ptr flags   ; set carry flag to what it was the
  136.     sahf                    ; last time we exited this routine
  137.  
  138.     mov cx,24               ; will do first loop 24 times
  139.     mov ax,ds
  140.     mov es,ax
  141.     mov si,offset(swbseed)+13*4; set up ds:si -> x[13]
  142.     mov di,offset(swbseed)     ; set up es:di -> x[0]
  143.  
  144. loop1:  ; On a 80386, the instructions `lodsd', `stosd' and `loop'
  145.     ; all change registers automatically as noted in parentheses
  146.  
  147.     lodsd                   ;    ax = x[i+13]     ( si = si+4 )
  148.     sbb eax,[di]            ;    ax = ax-x[i]-carry
  149.     stosd                   ;    y[i] = ax        ( di = di+4 )
  150.     loop loop1              ; loop for i=0..47    ( cx = cx-1 )
  151.  
  152.     mov cx,13               ; will do next loop 13 times
  153.     mov si,offset(swbseed)  ; set up ds:si -> x[0]
  154.                                 ; es:di is already set up
  155. loop2:
  156.     lodsd                   ;    ax = x[i-24]     ( si = si+4 )
  157.     sbb eax,[di]            ;    ax = ax-x[i]-carry
  158.     stosd                   ;    x[i] = ax        ( di = di+4 )
  159.     loop loop2              ; loop for i=48..73   ( cx = cx-1 )
  160.  
  161.     lahf
  162.     mov byte ptr flags,ah   ; save carry flag for next time
  163. ;
  164. ; XOR the elements of swbb with a congruential generator and put the
  165. ; result in swbnn.x, reset the counter and the pointer.
  166. ;
  167.  
  168.     mov  di,bx
  169.     mov  [bx-4],bx
  170.     mov  si,offset(swbseed)
  171.     mov  cx,N
  172. IFNDEF fast
  173.     mov  eax,congx
  174.     mov  ebx,69069
  175. loopc:  mul  ebx
  176.     mov  edx,eax
  177.     lodsd
  178.     xor  eax,edx
  179.     stosd
  180.     mov  eax,edx
  181.     loop loopc
  182.     mov  congx,eax
  183. ELSE
  184.     rep  movsd
  185. ENDIF
  186.     ExitFill
  187. SWBfill    endp
  188.  
  189. ; Random Number procedures ============================================
  190. ;
  191. CheckFill MACRO bits,bytes,count
  192.     local ok
  193.     dec    bits.c
  194.     jns    ok
  195.         mov   bx,offset(bits.x)
  196.         call  SWBfill
  197.         mov   bits.c,count-1
  198. ok:    mov    bx,bits.p
  199.     add    bits.p,bytes
  200. ENDM
  201.  
  202. i32bit proc
  203.     EnterProcedure
  204.     CheckFill swb32,4,N
  205.     mov    ax,[bx]
  206.     mov    dx,[bx+2]
  207.     DwordFn
  208.     ExitProcedure
  209. i32bit endp
  210.  
  211. i31bit proc
  212.     EnterProcedure
  213.     CheckFill  swb32,4,N
  214.     mov    ax,[bx]
  215.     mov    dx,[bx+2]
  216.     and    dh,7Fh
  217.     DwordFn
  218.     ExitProcedure
  219. i31bit endp
  220.  
  221. i16bit proc
  222.     EnterProcedure
  223.     CheckFill  swb16,2,2*N
  224.     mov    ax,[bx]
  225.     WordFn
  226.     ExitProcedure
  227. i16bit endp
  228.  
  229. i15bit proc
  230.     EnterProcedure
  231.     CheckFill  swb16,2,2*N
  232.     mov    ax,[bx]
  233.     and  ah,7Fh
  234.     WordFn
  235.     ExitProcedure
  236. i15bit endp
  237.  
  238. i8bit  proc
  239.     EnterProcedure
  240.     CheckFill  swb8,1,4*N
  241.     mov    al,[bx]
  242.     ByteFn
  243.     ExitProcedure
  244. i8bit  endp
  245.  
  246. i7bit proc
  247.     EnterProcedure
  248.     CheckFill  swb8,1,4*N
  249.     mov    al,[bx]
  250.     and al,7Fh
  251.     ByteFn
  252.     ExitProcedure
  253. i7bit  endp
  254.  
  255. i1bit  proc
  256.     EnterProcedure di
  257.     dec    swb1.c
  258.     jns    ok1
  259.         CheckFill  swb32,4,N    ; do an i32bit
  260.         mov      dx,[bx+2]
  261.         mov   bx,[bx]        ; with the answer in dx:bx
  262.         mov   swb1.c,31
  263.         mov   di,offset(swb1.x)
  264.         mov   swb1.p,di
  265.         mov   ax,ds
  266.         mov   es,ax
  267.  
  268.         mov   cx,16
  269. stosb1:        mov      al,1
  270.         and   al,bl
  271.         stosb
  272.         shr      bx,1
  273.         loop  stosb1
  274.  
  275.         mov      cl,16
  276. stosb2:        mov   al,1
  277.         and   al,dl
  278.         stosb
  279.         shr   dx,1
  280.         loop  stosb2
  281. ok1:    mov    bx,swb1.p
  282.     inc    swb1.p
  283.     mov    al,[bx]
  284.     ByteFn
  285.     ExitProcedure di
  286. i1bit  endp
  287.  
  288. uni proc
  289.     EnterProcedure
  290.     fild    neg31
  291.     CheckFill    swb32,4,N
  292.     and    byte ptr [bx+3],7Fh
  293.     jnz    oku
  294.         mov   eax,[bx]
  295.         mov   tmpdhi,eax
  296.         CheckFill  swb32,4,N
  297.         mov   eax,[bx]
  298.         mov   tmpdlo,eax
  299.         fild  tmpq
  300.         jmp   uxit
  301. oku:    fild    dword ptr [bx]
  302. uxit:    fscale
  303.     fstp    st(1)
  304.     RealFn
  305.     ExitProcedure
  306. uni endp
  307.  
  308. vni proc
  309.     EnterProcedure
  310.     fild    neg31
  311.     CheckFill    swb32,4,N
  312.     test    byte ptr [bx+3],0FFh
  313.     jnz    okv
  314.         mov   eax,[bx]
  315.         mov   tmpdhi,eax
  316.         CheckFill  swb32,4,N
  317.         mov   eax,[bx]
  318.         mov   tmpdlo,eax
  319.         fild  tmpq
  320.         jmp   vxit
  321. okv:    fild    dword ptr [bx]
  322. vxit:    fscale
  323.     fstp    st(1)
  324.     RealFn
  325.     ExitProcedure
  326. vni endp
  327.  
  328. duni proc
  329.     EnterProcedure
  330.     fild    neg63
  331.     CheckFill    swb32,4,N
  332.     dec    swb32.c
  333.     jns    okdu
  334.         mov   eax,swb32.x[4*N-4]
  335.         mov   tmpdlo,eax
  336.         mov   bx,offset(swb32.x)
  337.         call  SWBfill
  338.         mov   swb32.c,N-1
  339.         mov      eax,swb32.x
  340.         mov   tmpdhi,eax
  341.         and   byte ptr [tmpq+7],7Fh
  342.         fild  tmpq
  343.         jmp   duxit
  344. okdu:    and    byte ptr [bx+7],7Fh
  345.     fild    qword ptr [bx]
  346. duxit:    add    swb32.p,4
  347.     fscale
  348.     fstp    st(1)
  349.     DoubleFn
  350.     ExitProcedure
  351. duni endp
  352.  
  353. dvni proc
  354.     EnterProcedure
  355.     fild    neg63
  356.     CheckFill  swb32,4,N
  357.     dec    swb32.c
  358.     jns    okdv
  359.         mov   eax,swb32.x[4*N-4]
  360.         mov   tmpdlo,eax
  361.         mov   bx,offset(swb32.x)
  362.         call  SWBfill
  363.         mov   swb32.c,N-1
  364.         mov   eax,swb32.x
  365.         mov   tmpdhi,eax
  366.         fild  tmpq
  367.         jmp   dvxit
  368. okdv:    fild    qword ptr [bx]
  369. dvxit:    add    swb32.p,4
  370.     fscale
  371.     fstp    st(1)
  372.     DoubleFn
  373.     ExitProcedure
  374. dvni endp
  375.